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Abstract Time-distance helioseismology uses cross-covariances of wave motions on the 
solar surface to determine the travel times of wave packets moving from one surface loca¬ 
tion to another. We review the methodology to interpret travel-time measurements in terms 
of small, localized perturbations to a horizontally homogeneous reference solar model. Us¬ 
ing the first Born approximation, we derive and compute 3D travel-time sensitivity (Frechet) 
kernels for perturbations in sound-speed, density, pressure, and vector flows. While kernels 
for sound speed and flows had been computed previously, here we extend the calculation to 
kernels for density and pressure, hence providing a complete description of the effects of 
solar dynamics and structure on travel times. We treat three thermodynamic quantities as in¬ 
dependent and do not assume hydrostatic equilibrium. We present a convenient approach to 
computing damped Green’s functions using a normal-mode summation. The Green’s func¬ 
tion must be computed on a wavenumber grid that has sufficient resolution to resolve the 
longest lived modes. The typical kernel calculations used in this paper are computer in¬ 
tensive and require on the order of 600 CPU hours per kernel. Kernels are validated by 
computing the travel-time perturbation that results from horizontally-invariant perturbations 
using two independent approaches. At fixed sound-speed, the density and pressure kernels 
are approximately related through a negative multiplicative factor, therefore implying that 
perturbations in density and pressure are difficult to disentangle. Mean travel-times are not 
only sensitive to sound-speed, density and pressure perturbations, but also to flows, espe¬ 
cially vertical flows. Accurate sensitivity kernels are needed to interpret complex flow pat¬ 
terns such as convection. 
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1 The forward problem in time-distance helioseismology 


There are two general approaches to helioseismology. Global helioseismology consists of 
interpreting the frequencies of the normal modes of oscillations, whereas local helioseis- 
mology uses spatial-temporal correlations of the wave field as the input data (time-distance 
helioseismology, acoustic holography, direct modeling) or local power spectra of solar os¬ 
cillations ( ring-diagram analysis). For overviews of local heliose ismol ogy techniques, se e 
reviews bv IChristensen-Dalsgaard d2002l) . lGizon and Birchl (l2005h . and lGizon et al] J201dh . 

In this paper, we focus on the forward problem of time-distance helioseismology, i.e. 
how to compute the effects of subsurface heterogeneities on wave travel times. The problem 
of inferring subsurface heterogen eitie s from observed travel times (the inverse problem) was 
first formalised and solved by Kosovichevl ( 1199 6) and is also discussed in the reviews cited 
above. Throughout this paper we consider travel times measured between two points on 
the solar surface using a cross-covariance technique C Du vail et alJ 19931) . The purpose of the 
paper is to review and extend the computation of the sensitivity of wave travel times to small- 
amplitude changes to the local physical conditions of the solar interior. By a small-amplitude 
change to a reference solar model we mean a perturbation that causes a phase shift in the 
observed wave field that is much less than a radian. These sensitivity functions depend on the 
3D position in the solar interior (the coordinates of local scatterers) and are known as Frechet 
kernels, or simply kern els, in Earth seismology (Marquering et al .11 1 9981 . 1 1 999l : fDahlcn et alJ 
l200(l : IHung et al .12000! ). The solar problem differs from the earthquake seismology problem 
because solar oscillations are continuously excited by a random process, turbulent solar 
convection (granulation). 

The first kernel calculations in time-distance helioseismology employed the ray approx¬ 
imation, in whi ch the wavelength is assumed to be small compa r ed to the other length scales 
in the problem d Kosovichevl 1 996l : lKosovichcv and Duvalll 1 997h . [Bi rch and K osovichevl d2000l) 
used the first Born approximation (single-scattering approximation) to compute 3D ker¬ 
nel^ for sound-speed perturbations using the single-source approximation. Subsequently, 
iGizon and Birchl c2002l) introduced a general framework for Born kernel calculations includ¬ 
ing a treatment of wave excitation by near-surface convection, this framework is the basis for 
the calculations presented here. In particular. IGizon and Birch (2002 ) showed that spatially 
distributed wave sources alter the sensitivity kernels: scattering far away from the observa¬ 
tion points can be important. They also showed that the details of the wave power spectrum, 
wave attenuation, and the details of the measurement procedure must all be included in 
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the kernel computations. They proposed a measurement procedure that enables a consis¬ 
tent computation of kernels by relating small-amplitude changes in the cross-covariance 
to travel-time shifts. Usin g this approach ke rnels have been computed for source amph 


2004), and vector flows JOizon et al. 200(1 


tude and wave attenuation (Gizon and Birch 2002), sound-speed perturbations (Birch et al 


Jackiewicz et alJ1200a . 12007 ; Birch and Gizon 


2007; Birch et al. 2 011). We note that th e same approach was used to compute kernels for 


ring-diagram analysis (Birch et al.ll2007l) . All of these calculations assume that the helio¬ 
seismic observable is the vertical wave velocity at the solar surface. Away from disk center, 
the horizontal velocity contributes to the Doppler velocity and short wavelengths are less 
visible towards the limb lijackiewicz et al.lF2006 l. 1 2007 ). All of these calculations employed 
a plane-parallel approximation for the background solar model (used for the study of small 
patches on the Sun). 

More recent approaches based on the adjoint method are particularly useful for solving 
problems with 3D heterogeneous background models lHanasoge et al.ll201 ll) . This method 
has been applied to a background model containing a magnetic field and 3D kernels_for 
sound-speed, density, 3D flows, and magnetic fields were computed (iHanasocc et al.l2012l) . 
Mathematic al expressions for density kernels, as w ell as sound-speed kernels, have also been 
presented bv ISchlottmann and Kosovichevi ( j20T~fl) using “functional analytic methods”. The 
contribution of the lowest-order non-zero perturbat ions of magnetic field, to a background 
containing no magnetic field, is discussed in lGizon and Thompsonl d2007|). 

Ot her approaches than the Bom approximation exist. For example, fjensen and Piipersl 
120031) applied the Rytov approximation, which is closely related to the Born approximation, 
to compute kernels for sound-speed and flows in the single^source approximation. 

In this paper, we work within the framework of iGizon and B ire hi |2002l) and extend the 
computation of travel-time kernels to treat density and pressure perturbations, in addition 
to sound-speed and flow perturbations. For the sake of computational simplicity, we work 
in Cartesian geometry starting from a ID static and non-magnetic reference solar model. 
Furthermore, we do not assume a priori hydrostatic equilibrium in the perturbed model to 
allow for contributions to the force balance from flows and magnetic fields and thus we 
retain three thermodynamic variables (sound speed, density, and pressure in the perturbed 
model). 


2 Equivalent-source description of wave interaction 

Linear adiabatic stellar oscillations (e.g. iLvnden-Bell and Ostrikcrlil 9671) are typically gov¬ 
erned by equations of the form Jf § (r, t) = S(r, /), where <5? is a linear differential operator, 
£(r,f) is the vector-valued wave displacement, and S(r,t) is the vector-valued source func¬ 
tion responsible for wave excitation. Scattering events are uniquely identified by both time t 
and a three dimensional (3D) position vector r. 

In the first-order Born approximation (e.g. lGizon and Birchl2002l) . the standard coupled 
system of linear differential equations is 

2%So = So, (1) 

Jf 0 S% = -SJf^ 0 +SS. (2) 

Equation 0 describes the propagation of small amplitude waves through a background 
solar model. The operator is the reference linear wave operator, E, (] the displacement 
field in the reference model and So is the reference source field. Equation {2} governs the 
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perturbed displacement field, The first-order linear operator 8,'Z J encapsulates pertur¬ 
bations to the solar thermodynamic structure and the introduction of flows, and 5S is the 
perturbation to the source field. In this paper, we ignore 8S (see e.g. iGizon and Birchl[2002l . 
for a discussion). 

For a reference solar model with vanishing flows and magnetic fields, and also using the 
Cowling approximation, the reference linear operator governing the reference displacement 
field is (see e.g. ICameron et alj|2006l) 

-^olo := Po {d t +r) 2 § 0 — V (cqPo V- | 0 + lo' ^ p o) +go V • (po4o) ■ (3) 

The background solar model variables, which depend on the distance from the center of the 
Sun, have their usual symbols: Co is the sound-speed, po is the density, Po is the pressure, 
and go is the acceleration due to gravity. The partial derivative with subscript t denotes 
differentiation with respect to time and the vector V is the 3 dimensional spatial gradient 
operator. The phenomenological operator, F, is responsible for damping the oscillations. 
The random source function So is assumed to be statistically homogeneous in the horizontal 
direction and statistically stationary in time and represents the driving of oscillations by 
near-surface turbulent convection. 

The first-order linear operator due to first-order perturbations in sound-speed squared 

(8c 2 ), density (8p), pressure (8P) and 3D mass flows (v) is, 

5-S^o := (8c 2 poV ■ %o) +2p 0 (v-V)<9,£ 0 — V(| 0 • VSP) 

+8p (8,+D 2 4 0 -V (c 2 8p V ■ | 0 ) + g 0 V • (8p | 0 ) + 8g V • (p 0 £„)• (4) 

The first-order changes in gravity 5g are related to first-order changes in density and are 
expected to be very small near the solar surface. 

We use sound-speed, density, and pressure as our three independent thermodynamic 
variables. It is possible to make other choices; e.g. the adiabatic exponent y may be intro¬ 
duced at the expense of the others by using the relationship p c 2 = Py. We do not assume 
hydrostatic equilibrium for the thermodynamic perturbations. If hydrostatic equilibrium was 
assumed, then it would be possible to eliminate the gradient of the pressure perturbation (in 
Eq. PI ) with a contribution from the density perturbation (or vice-versa), and consequently, 
reduce the number of independent perturbation vari ables to two (e.g. sound-speed and den¬ 
sity) as is done in global helioseismology (e.g. [Christensen-Dalsgaardj[2002l) . 


3 Computing travel-time sensitivity kernels for structure and flows 

In time-distance helioseismology wave travel times (t) are measured from the temporal 
cross-covariance of the observed signal, 0(x,t): 

1 f T 

C(f;xi,x 2 ) = - / )0(x 2 ,f +t), (5) 

T Jo 

where t is the time lag, xi and x 2 are two spatial locations on the solar surface and T is the 
duration of the time series. In the temporal Fourier domain it becomes 

C(ft>;xi,x 2 ) = h m 0*(xi ,co)(t>{x 2 ,co), (6) 

where co is the angular frequency and the superscript * deno tes complex_conjugate (through¬ 
out this paper we use the Fourier conventions from lGizon and Birchll200j) . C(©;xi,x 2 ) is 
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the temporal Fourier transform of C(f;xi,X 2 ) and (j)(x,co ) is the temporal Fourier trans¬ 
form of (j>(x,t). The sampling in frequency space is inversely proportional to the duration 
ha — /T. 

The observable <p is often a line-of-sight Doppler velocity; it is connected to the oscil¬ 
lations of the atmosphere caused b y sei smic waves. There exists a non-trivial relationship 
between | and (j) (see e.g. lNagashima et aJ]l2Q14l ) which can be expressed in general as 

<t>{x,co) = 0{|(r,©)}, (7) 


where O is an operator that acts on the 3D displacement and maps it onto the CCD detector. 
The operator includes the Point Spread Function (PSF) of the instrument plus any filters that 
were applied to the data. 

Then using the Born approximation (Section [2} the first-order change to the cross¬ 
covariance is given by 

5 C(©;x 2 ,xi) = —( 27 r) 4 £ f r/r£^(x 2 ,©,r) e r X(r, ©,xi) + (1 o 2 )*, ( 8 ) 
a 7© j 

= Yj[ *^ a (r,ffl;x2,xi)59 a (r) + (l© 2 )*, ( 9 ) 

a J® 

where X(r,a>,xi) := £ [ dr s [ Jr'(f'*(xi,co,r s )G''(r,co,r')M // ,(r J ,r',®), (10) 

nf 

and M u i{r s ,r' s , co) := E[S*{r s , a))5 , ; /(r',co)]. (11) 


The Green’s tensor G' is defined by 


2z?oG'(r,r-r s ,r s ) = <5(r-r s ) 5(r-t s )e ; , 


( 12 ) 


and for convenience we also defined 


«f' = 0{ G'}. (13) 

The operator 5^.“^ is the component of the perturbed wave operator (5Jzf) that is due to the 
change in the physical quantity 5q a (r). The quantity M,-,-/ is the source covariance tensor. 
The symbol (1 -B- 2)* is a short notation to mean switch xi and x 2 and take the complex 
conjugate of all other terms in the expression. The scattering physics is included in the 
functions % that give the linear sensitivity of the expectation value of the cross-covariance 
to changes in the solar model. In order to arrive at Equation ((Dj), it is often required to use 
integration-by-parts to separate the perturbation variables from spatial derivatives, and thus 
in those cases it is been assumed that relevant quantities vanish on the boundaries. 

Using Equation i[9j. '6" a can be expressed for each of the perturbations 


^ e 2(r;x 2 ,xi) 


%(r;x 2 ,xi) 


— (27T) 4 p 0 (r)c5(r) 


X>.^(x 2 ,r) 


V r ■ X(r, Xl ) + (1 h2)*, 


— (27r) 4 p 0 (r) [(co + iy) 2 £(f'(x 2 ,r)X / (r,xi) 
j 


+co(r) V r -X(r,x 1 )^a, 7 ^-'(x 2 .r) 
j 


(14) 


-X(r, Xl ) • V r £cf'(x 2 ,r)g ; -(r)j 
j 
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-(2 K) A 8 gi (r) ( r ) Sf (x 2 ,r) V r ■ [p 0 X(r, Xl )] + (1 o 2)*, 


^p(r;x 2 ,xi) 

^v,(r;x 2 ,xi) 


(2^) 4 P 0 (r)V. 


X(r,xi)£a, v ^(x 2 ,r) 


+ (1h2)‘ 


(2^) 4 i2®po(r)^- / (x 2 ,r)<9 r X 7 (r,xi) + (l o 2)* 


(15) 

(16) 
(17) 


where i is the imaginary number. Note that the first-order perturbation to gravity can be 
expressed in terms of first-order density according to 


V • §g = -47rG5p 


(18) 


where G is the gravitational constant. 

The integrals in Equations ©-Go} are 3D and over the volume of the Sun. The 3D flow 
vector (v) has been decomposed into its three components denoted v, . The sensitivity kernel 
K a for a given physical quantity q a is computed with all other physical quantities held fixed. 
Each kernel K a (r) due to a point scatterer at r in the interior is computed according to 


X„(r;xi,x 2 ) = 47rRe / VT*(®;xi,x 2 )^„(r,©;xi,x 2 )£/®, (19) 

Jo 

where Re takes the real part of the following expression and the function W gives the linear 
dependence of the travel-time shifts on the cross-covariance^The functions W incorporate 
all the details of the travel-time measurement procedureJ Gizon_and_Bircl]l2002h . One-way 
travel-times are sensitive to all types of perturbations. iDuvall et al. 1 19971 ) introduced travel¬ 
time differences cJTdjff = 5t(xi,x 2 ) — 5t(x 2 ,xi) to partially isolate the contribution from 
flows and mean travel times <>T mean = [5t(xi,x 2 ) + 5t(x 2 ,xi)]/2 to measure the other con¬ 
tributions. The corresponding W-functions are 


W d iff(©;x 1 ,x 2 ) = 1V(©;xi,x 2 ) -W(©;x 2 ,xi), (20) 

W mean (©;xi,X 2 ) = [W(©;xi,x 2 )-HT(©;x 2 ,xi)]/2. (21) 


The forward problem consists of modeling the chain of perturbations 8% —> 8(j) —> 5C —» 
8 T caused by changes in the structure and dynamics of the solar interior. Travel-time shifts 
(5 t) due to first-order Born perturbations in sound-speed squared, density, pressure and 3D 
vector flows are given by 


§t(xi,x 2 ) = / 
Ja 


= / dr 

0 


8 c ^ 8p 

K c 2 {r;\ u \ 2 )— (r)+Xp(r;xi,x 2 ) —(r) 
c o Po 

8P 3 

+Kp( r;x! ,x 2 ) — (r) + £ K Vi { r;x! ,x 2 ) v,-(r) 
rt) , = 1 


( 22 ) 


4 Plane parallel example 

Here we compute example kernels assuming a plane-parallel model which is appropriate for 
studying small areas on the Sun. We assume the observable is given by 

0(k, co) = —ia>F(k, co)£ ^ (k,z obs , (o), (23) 

where k is a horizontal wave vector, t is the line-of-sight unit vector that we assume points 
in the vertical direction (£ = z). We used an observation height of z 0 b s = 300 km. 
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The function F is a combination of both the instrumental PSF and filters. In this paper, 
we focus on modeling observations of waves with the same radial order. The solar power 
spectrum shows ridges of enhanced power. Each of these ridges corresponds to waves with 
a particular radial order. To isolate a single radial order, the power s pectr um is filtered using 
a ridge filter that suppress all waves at other radial orders (e.g. see lJackiewicz et ai]l2008h . 

In this paper we consider both an f-mode ridge filter and p \-mode ridge filter. 

We assume spatially uncorrelated sources, such that the example source covariance ten¬ 
sor is 

Mui(r s ,r' s ,co) = {2n) 2 8i,8^ z 8{x s -x' s )m((o)d Zs 8{z s -z S rc)d z ' s 8 (i-Zac), (24) 

where r s = ( x s ,z s ) and r' = (x',z'). This is the source covariance described bv lBirch et alJ 
d2004l l. The depth dependence is a vertical derivative of a vertical momentum source. This 
form of the source covariance gives a good match to the observed power spectrum. The fre¬ 
quency dependence here is Gaussian as in lBirch et alJ (|2004) h owever, there is an additional 
factor of l/co 2 since we are solving for displacement whereas [Birch et aD ( 20041) was solv¬ 
ing for velocity. Here, the source auto-correlation function is in(co) := ft)~ 2 exp[—(ft>r sre ) 2 ], 
with a correlation time T SIC = 48 s. We note that M is defined to within an overall scale 
factor, which changes the amplitude of the ^ but does not change the kernels K q since 
the W-functions are inversely proportional to this scale factor. We used a source height of 
£src = 100 km dBirch et ai]|2004r) . 

4.1 Convenient mode-summation formula for the Green’s tensor 

The solutions to Equations (03 and (2j can be written in terms of Green’s functions ( jGizon and Birchl 
I 2 OO 2 I). We use normal-mode summation to compute Green’s functions, however, the ap¬ 
proach presented here is different to what was implemented in [Birch et ah] <2004). Here 
we write the Green’s functions for the damped problem in terms of normal modes of the 
undamped problem, and consequently, the new approach is less mathematically and com¬ 
putationally demanding. The equation governing the reference solar oscillations can be ex¬ 
pressed as 

-S?| 0 = Po [W +r) 2 + Jf] | 0 = So, (25) 

where can be identified from Section [2] Normal-mode eigensolutions (undamped) thus 
must satisfy 

jr£ n {r,k) = a&(k)S n {r,k), (26) 

where % " (r, k) is the eigenfunction with radial order n and horizontal wave vector k, co n ( k ) is 
the associated (real) eigenfrequency, and k := | |k| | is the wavenumber. Notice that because of 
the assumed horizontal isotropy of the reference solar model the normal-mode frequencies 
are independent of the direction of k. The displacement eigenfunctions may be decomposed 
into vertical and horizontal components according to 

I" (r,k) = [z^(z,k)+ik^(z,k)\e' kx , (27) 

where x is a horizontal position vector and z denotes height. Because the only restoring 
forces are pressure and buoyancy there is no motion in the direction k x z. We choose to 
normalise the eigenfunctions according to 

/ Po£"*(r,k) • ^' !, (r,k') rf 3 r=5„ y 5(k-k'), 


(28) 
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where 8 nn i is a Kronecker delta, 5(k —k') is a Dirac delta function. The integral is three 
dimensional and over all space. For the horizontally homogeneous medium under study, 
we define the Green’s vector G'(x —s ,z,t — t s ,z s ) as the displacement vector at position 
r = (x,z) and time t that results from a delta function source at position r s = (s,z s ) and time 
t s pointing in the e,- direction, where e, is one of the unit vectors e x , e v or e z . If we take the 
Fourier transform in time, then we get 

po[-(co + iy) 2 + Jf’]G‘(x-s,z,co,z s ) = ^-5(r-r s )e,-. (29) 

171 

The function y damps the waves so their lifetimes are finite. Typically it is chosen such 
that the subsequent line-widths of the model reproduce the line-widths observed in the Sun. 

We now expand the G' on a basis of undamped normal-mode solutions according to 

G'(x-s, z, w, z s ) = Y< [ ^ 2 kc^(ft»,r s ,k)|"(r,k). (30) 

n J 


where cJ,(o>,rs,k) are coefficients to be determined. The sum is over all radial orders n 
and the integral is two dimensional over all horizontal wave numbers. After taking the dot 
product of Equation ( 1291) with *(r,k') and integrating over r, we find 


[- ((o + i /„/ (k) ) 2 + co „ 2 , (k)} cj,, ( co, r s , k) = — <§/’ * (r s , k) , 


(31) 


where y(k) £"(r,k) = y,(k)%"( r.k). Finally, after some additional manipulation we have 


G'(k,z,©,z s ) 


J_ C (^;k) ff*(z s ,k) 

2^r-[®+i %(*)] 2 +«£(*)’ 


(32) 


where we introduce the convenient notation for the eigenfunction 

C"(z,k) :=z^(z,k)+ik^{z,k). (33) 

Equations d32b and d33b give a convenient set of expressions for computing the Green’s 
functions for the damped problem in terms of the modes of the undamped problem. They 
are, in practice, easier to use than the expressions of lBirch et al.l (120041 ). which involve the 
modes of the damped problem; the modes of the undamped problem are easier to compute 
than the modes of the damped problem. In addition, the damping model (i.e., the choice of 
the y„) can be updated independently of the computation of the normal modes. 

Figure [I] shows how the complex Green’s function presented here compares with two 
computations from iBirch et ahl (2004) versus acoustic depth, at fixed co and k. iBirch et all 
120041 ) used a slightly different definition of the Green’s function: (d, + F) 2 is replace d by 
d; 2 + 2 F d , in the wave operator ,!// and F = F(z). The problem was solved by IBirch et al.l 
( 2004 ) by summation over the modes of the damped problem (unlike in the present pa¬ 
per) and also by direct numerical solution of the differential equations. The source depth is 
z s = 3.7 Mm and the horizontal wavenumber is k = 1 radMm~*. For the 2004 paper, the 
cyclic frequency is co/2k = 3.92 mHz, a frequency just above the n = 2 resonance. The 
normal mode eigenfrequencies from the model used in this paper are slightly different from 
those in iBirch et al.l (1 20041) : on the order of a few uHz. Thus in order to make a good com¬ 
parison the new Green’s function was evaluated a few /iHz higher at co/ 2k = 3.92265 mHz, 
in order to be a similar distance from the n = 2 mode. We find that the two mode-summation 
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solutions (both used radial orders 0 < n < 14) give similar answers (dashed curves in Fig¬ 
ure [T]l: The imaginary parts of the Green’s functions match very well, the real parts less 
so. The difference in the real parts_is due to the non-vanishing imaginary component of 
the damped eigenfunctions from lBirch et alj (20041) . The real part of the numerical Green’s 
function (solid curves in Figure Q} has a discontinuity at the source depth that cannot be 
reproduced by a sum over only 15 modes. Whether the above differences are important for 
the interpretation of helioseismic travel times deserves to be studied. 


x 10' 5 



Fig. 1 Horizontal component of the Green’s function G z fro m this study (blue dashed curves) compared with 
normal-mode summation Green’s function from lBirch et alj d20041) (dashed black curves) and the numerical 
solution from iBirch et ahl 120041) (solid black curves). There are two curves for each Green’s function to 
show the real and imaginary parts. The source depth is Zs — 3.7 Mm and the horizontal wavenumber is 
k = 1 rad Mm -1 . The mode summation formulae uses modes from n — 0 up to n — 14. The vertical black 
dashed line shows the observation height and the vertical solid black line shows the source depth. 


4.2 Background Solar Model and Zero Order Wavefield 

We use model S dChristensen-Dalsgaard et alii 19961) for the background model of the Sun, 
where sound-speed, density, pressure, and gravity are functions of height only. There are no 
flows or magnetic fields in this background model. 

The reference displacement field for th e undamped problem (i .e. £ n ) is computed by 
solving Equation {TJ. We follow the work of lGizon and Birchl |2002) and I Birch et a l. (2004) 
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and use a plane-parallel (Cartesian) approximation. The Cowling approximation^ also en¬ 
forced. This reduces to an eigenvalue problem similar to that encountered in iBirch et alJ 
( 12004 ). however here we solve for rather than df% 0 . It was particularly important to en¬ 
sure the eigenfunction solutions were well behaved under double numerical differentiation 
over height, since some of the features of the pressure kernels are sensitive to small errors. 
We use the same source covariance as in lBirch et al.l (2004). 

The undamped reference displacement field is used to compute the attenuated Green’s 
function through Equation d32b . 


4.3 Numerical resolution in Fourier space and spatial computational domain 

The kernels were computed using the Model-S non-uniform height grid that is provided 
as part of the Aarhus adiabatic pulsation package. We selected to use a subset of the height 
domain from z = 0.4959 Mm to z — —20.9066 Mm. We use a uniform horizontal grid (Ax = 
Ay = 96 km) spanning 96 Mm in each direction. The size of the computational box was 
selected to ensure the kernels essentially vanish on the horizontal boundaries and also the 
bottom boundary. In particular, the amplitudes of the hyperbolic features of the kernels are 
notorious for decaying slowly with distance. We conducted resolution tests to confirm the 
example kernel s prese nted in this paper did not change significantly with higher resolution. 

Fo llowing lBirch et al.l ll20041. we used the line-width measurements from Korz ennik et al.l 
J2004I) as a basis for the empirical damping model. These line widths set the numerical 
resolution in Fourier space that is required in the calculation of the kernels. For the p\- 
mode examples shown in this paper, we used frequency and wavenumber resolutions of 
Av = 3.8 liHz and Ak = 1.59 x 10~ 4 rad Mm -1 respectively. The frequency domain ranged 
from / = 1.5 mHz to / = 5.3 mHz. The wavenumber domain was from k = 0.11 rad Mm~* 
to k = 2.5 rad Mm -1 . We computed higher resolution tests to confirm that the resolution 
used here was sufficient. 

To obtain a good representation of p\ modes (including the line asymmetry), normal 
modes from = 0 to n = 14 were used in the summation for the Green’s functions. A filter 
was then applied to isolate the p\ ridge. We used a Modulation Transfer Function (MTF) of 
the form e~^ k where j6 = 1.15 Mm. 


4.4 Example kernels 

The code used to compute travel-time sensitivity kernels for this paper was writte n in MAT- 
FAB and is a modified and extended vers ion of the code used by IBirch et all d2004h to 
compute mean sound-speed kernels and bv IBirch and Gizonl (1 20071) to compute difference 
velocity kernels. The main modifications are the computation of density and pressure kernels 
and the new method to compute Green’s functions (1 14.11 ). For simplicity, when computing 
the density kernel we neglected the perturbation the gravitational field (5g) and yin Equa¬ 
tion HU: both these effects are small. 

Figure[2]shows horizontal and vertical slices through kernels for the changes in the mean 
travel-time due to fractional perturbations to the squared sound-speed, density, and pressure 
(for pi-modes) and a separation of A = X 2 — xi = 16.0347 Mmx between the observation 
points. All the kernels are symmetric in both horizontal coordinates and show the same 
general features as the sound-speed kernels shown bv lBirch et alJ (200 41. The kernels have 
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complicated features and are largest in amplitude near the two observation points. The ker¬ 
nels extend roughly over the first few Mm below the photosphere; this is the region where 
the pi kinetic energy density is significant for frequencies 1.5-5.3 mHz. The total integral of 
the sound-speed kernel is negative (increased sound-speed leads to a reduction in <5T mn ). The 
total integral of the density kernel is also negative implying a uniform increase in density 
would also lead to a reduction in (5 r nln . Conversely, the total integral of the pressure kernel 
is positive, which implies a uniform increase in pressure would lead to an increase in §T mn . 
The pressure and density kernels are strongly anti-correlated. 



Fig. 2 Mean point-to-point travel-time sensitivity kernels for relative changes in sound-speed squared (top 
panels), density (middle panels), and pressure (bottom panels). Kernels have units of s Mm -3 . Normal modes 
n = 0 to n = 14 were used to compute the Green’s function and a filter was applied to isolate the p \-mode 
ridge. The left column displays horizontal cuts at a height of z = —0.3 Mm. The right column displays vertical 
slices along the line y — 0 Mm. All colour scales are symmetric about zero (green) with red positive and blue 
negative. The crosses in the left panels represent the two points Xi and X 2 ;the distance between the two points 
is 16.0347 Mm. 


Figure[3]shows slices through the kernels for the sensitivity of the travel-time differences 
to changes in sound speed, density, and pressure, for the same geometry and filter as the 
kernels shown in Figure [2] These kernels, unlike in Figure [2] are anti-symmetric in x, but 
still symmetric in y. This is a consequence of the definition of the travel-time difference (the 
difference between the two one-way travel times). The anti-symmetry in x is important as 
it implies that the total integral of these kernels vanishes, i.e., a spatially uniform change in 
solar structure does not produce a change in the travel-time differences. 

Figure[4]shows slices through the kernels for the sensitivity of the travel-time differences 
to changes in v x , v y , and v z , for the same geometry and filter as the kernels shown in Figure[2] 
The kernel for v z is anti-symmetric in x and symmetric in y, whereas the kernel for v y is 
anti-symmetric in both x and y. This implies that uniform perturbations in v z and v y will 
not effect travel-time differences. The difference kernel for v x is symmetric in both x and 
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Fig. 3 Difference point-to-point travel-time sensitivity kernels for relative changes in sound-speed squared 
(top panels), density (middle panels), and pressure (bottom panels). Kernels have units of s Mm -5 . Normal 
modes n = 0 to n = 14 were used to compute the Green’s function and a filter was applied to isolate the p\ - 
mode ridge. The left column displays horizontal cuts at a height of z = —0.3 Mm. The right column displays 
vertical slices along the line y = 0 Mm. The crosses in the left panels represent the two points xi and X 2 ; the 
distance between the two points is 16.0347 Mm. 


y and consequ ently has a non-zero to tal i ntegral. 3D differen ce kernels for the f-mode was 
presented in iBirch a nd Gizonl 1 2007) and Birch et alj 1 2 01 lh. The present calculations are 
fully consistent with Birch and Gizon ( 2007) and Birc h et alJ d201 lh . 


Figure [5] shows slices through the kernels for the sensitivity of mean travel-times to 
changes in v x , v y , and v z , for the same geometry and filter as the kernels shown in Fig¬ 
ure [2] The symmetries imply that uniform perturbations in both v x and v y would not affect 
mean travel time s. Equ ivalent 2D kernels for v x and v y , for the f-mode, were shown in 
Ijackiewicz et alj 1200 7'). We show here that the mean kernel for v z is symmetric in both x 
and y and has a significant non-zero total integral. 


Figure [6] demonstrates how the total integrals of the kernels converge as a function of 
distance from the observation points. The total integrals for the structure and horizontal flow 
kernels have approximately converged over the 96 Mm 2 domain, however, the mean kernel 
for vertical velocity has not yet converged. 

Table |H shows the mean and difference total integrals for p\ modes. This gives an indi¬ 
cation of how travel-time shifts respond to uniform perturbations. The largest contribution 
to the mean travel-times is from a uniform change in relative sound-speed. The contribu¬ 
tions to mean travel-times from relative density and pressure perturbations are almost equal 
but opposite in sign. For the flows, the only non-zero contribution to the mean travel-times 
would be a uniform vertical flow, the horizontal components vanish due to their symmetries. 
The only contribution to the difference travel times due to uniform perturbations is due to 
flows along the two observation points. 
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Fig. 4 Difference point-to-point travel-time sensitivity kernels for v x (top panels), v y (middle panels), and 
vz (bottom panels). Kernels have units of s Mm -3 / (km/s). Normal modes n = 0 to n = 14 were used to 
compute the Green’s function and a filter was applied to isolate the /?i-mode ridge. The left column displays 
horizontal cuts at a height of z — —0.3 Mm. The right column displays vertical slices along the line y — 0 
Mm. All colour scales are symmetric about zero (green) with red positive and blue negative. The crosses in 
the left panels represent the two points xi and X 2 ; the distance between the two points is 16.0347 Mm. 



JK c2 dr 

JKpdr 

jK P dr 

jK Vx d r 

jK Vy d r 

jK v .dr 

mean 

-460 s 

-170 s 

168 s 

0 s/(km/s) 

0 s/(km/s) 

54 s/(km/s) 

diff 

0s 

0 s 

0 s 

-152 s/(km/s) 

0 s/(km/s) 

0 s/(km/s) 


Table 1 Total integrals of the mean and difference kernels for sound-speed, density, pressure and flows for 
pi modes. 


5 A consistency check 


In this section, we compute the change in the cross-covariance (e.g. see Equation AC, 
caused by a simple sound-speed perturbation of the form 


Ac 2 

—=- = e exp 


-(z-zq ) 2 

2a 2 


where e = 0.0001, zo = —1.5 Mm, and a = (1 Mm)/(81n2) 1 / 2 . 
Because e is very small, AC can be approximated by 


r §c 2 

8C(co)= dxdz% 2 (x,z,co;x i,x 2 ) —(z) 

J Cr. 


(34) 


(35) 


where the two points are separated by a distance of 16 Mm. The function SC(w) is real. 
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Fig. 5 Mean point-to-point travel-time sensitivity kernels for v x (top panels), v y (middle panels), and vz 
(bottom panels). Kernels have units of s Mm~ 3 / (km/s). Normal modes n — 0 Lo n — 14 were used to compute 
the Green’s function and a filter was applied to isolate the p\ -mode ridge. The left column displays horizontal 
cuts at a height of z = —0.3 Mm. The right column displays vertical slices along the line y = 0 Mm. All colour 
scales are symmetric about zero (green) with red positive and blue negative. The crosses in the left panels 
represent the two points Xj and x^: the distance between the two points is 16.0347 Mm. 


Because Ac 2 depends on height only, we can also compute AC using 

AC(co) = C(co) —Cq(co). (36) 

where C(co) is the cross-covariance computed using the eigenmodes of the solar model 
with sound speed co + ^c- To compute the new eigenmodes we use the same solver as 
in Section 2.2. Here we used a different damping model where each mode was damped 
according to y„(k) = y|fflo(k)/©| 4 ' 4 , where m/(27t) = 3 mHz and f/(2.7t) = 100 Hz. In 
addition, for wave numbers less than 0.8 Mm we set the damping to remain constant. The 
resolution in wavenumber was the same as earlier, however the resolution in frequency was 
changed to 19.1 mHz. In Figure [7) we find that there is a very good agreement between 
AC(co) and 8C(w), which is expected in the limit e —> 0. These results give us confidence 
in the reliability of the computations. Furthermore, it appears that the numerical resolutions 
and the extent of the horizontal domain are sufficient. 


6 Are there simple relationships between the kernels? 

In order to produce reliable time-distance inversions of the solar interior it is helpful if the 
kernels have different structures and features. Figures[8]and[9]display cuts through the mean 
and difference kernels at z = —0.3 Mm along the x and y axes. We have plotted the negative 
values of the pressure kernel (green) in order to better illustrate the strong anti-correlation 
with the kernel for density (red). It can evident that kernels for sound-speed squared (blue), 
density (red) and negative pressure (green) have similar features. 
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Fig. 6 Total 3D integrals of the kernels as a function of A\ +A 2 — A (Mm). The flow kernels were scaled 
by a uniform flow of 1 km/s such that all the total integrals have units of seconds. The left panel is for mean 
kernels and the right panel for difference kernels. The integrals for sound-speed squared is blue, density is 
red, pressure is green, v x is cyan, v y is magenta, and v z is black. The total integrals for the mean kernels 
for sound-speed squared, density and pressure are close to converged at large distances from the two points 
however, for the vertical velocity a larger domain is necesssary to ensure convergence. The total integrals of 
the mean kernels for horizontal flows vanish. The only non-zero total integral for the difference kernels is the 
one for v x (cyan) and its total integral is close to converged with the domain of 96 Mm. 


Travel-time shifts due to sound-speed perturbations are directly related to the divergence 
of the eigenfunctions. Therefore, travel-times shifts for f-modes are not sensitive to sound- 
speed perturbations since the divergences of the f-mode eigenfunctions vanish. However, the 
power spectrum in the neighbourhood of the f-mode ridge has a small contribution from the 
wings of the p modes and thus f-mode mean travel-time measurements are weakly sensitive 
to sound-speed perturbations. For example, we find that the total integral of K c i for the f- 
mode ridge at A = 16.0347 Mm is approximately 3% of the total integral of K c 2 for the p\ 
ridge. 


7 Specific case of hydrostratic equilibrium 


In the fully general non-hydrodynamic case, where we have three independent thermody¬ 
namic perturbation variables, the travel-time shifts (ignoring flows) are given by 


5 1 = / dr 


L 


K 


(r\ 


Sc 2 


(37) 


We have now explicitly expressed which quantities are held fixed. The kernel K c i pP is 
due to perturbations in sound-speed at fixed density and pressure, similarly K pc i P is the 
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Fig. 7 Change in the cross-covariance (at a fixed point-to-point distance of A = 16.0347 Mm) for the case 
of the Gaussian perturbation to the sound-speed of Model S. The red curve is SC and the blue curve AC. The 
results show very good agreement. 


density kernel at constant sound-speed and pressure and finally K Pc i p is the pressure kernel 
at constant sound-speed and density. 

In global helioseismology, it is typically assumed that the density and pressure pertur¬ 
bations are in hydrostatic equilibrium (e.g 'Gough and Thompson!! 19911 ) 

V<5P = g 0 Sp, (38) 


where we assume the perturbation to the gravitational field is negligible for the sake of 
simplicity. 

Equation ( 138b may be manipulated to express the pressure perturbation entirely in terms 
of a density perturbation (or vice-versa) and thus the number of independent thermodynamic 
variables is reduced from three to two. 

The horizontal components of Equation d38b implies that that the pressure perturbation 
is at most a function of height, i.e. 8P = 8P(z). The vertical component implies that the 
density perturbation is also at most a function of height i.e 8p = 8p(z), and hence 


d8P 

—— = —go op. 


(39) 


It is possible to implement the assumption of hydrostatic equilibrium, for example by 
eliminating pressure in terms of density. To achieve this it is necessary to manipulate the 








Interpretation of Helioseismic Traveltimes 


17 




Fig. 8 Cuts through the mean point-to-point travel-time sensitivity kernels for sound-speed squared (blue), 
density (red), negative pressure (green), v x (cyan), v y (magenta), and vz (black). The sound-speed, density 
and pressure kernels have units of s Mm -3 , whereas the velocity kernels have units Mm -3 / (km/s). Normal 
modes n = 0 to n = 14 were used to compute the Green’s function and a filter was applied to isolate the 
pi-mode ridge. The cuts are at taken at a height of z — —0.3 Mm. The left panel is a cut along the x axis, 
whereas the right panel is a cut along the y axis. The distance between the two points is 16.0347 Mm. 


last term on the right-hand-side of Equation d37b as a gradient of the pressure. So since the 
perturbations are all functions of height only, we require 

dA (r) K p^ P (. r ) 

* Po(z) ( ’ 


with boundary condition A(x,zt) = 0. This could be solved for each x. Provided a solution 
exists, after integration by parts once, then the travel time shifts become 


5 T = 


Sk z. 

J 3 dr K c 2 , P p( r )-pr(z) + ( K P ,M r ) +Po(z)«o(zM(r)) 


sp, 

Po 


(41) 


However, since we know we have two independent thermodynamic quantities then it is also 
possible to express the travel time shifts in terms of a different set of kernels according to 



(42) 
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Fig. 9 Cuts through the difference point-to-point travel-time sensitivity kernels for sound-speed squared 
(blue), density (red), negative pressure (green), v x (cyan), v y (magenta), and vz (black). The sound-speed, 
density and pressure kernels have units of s Mm -3 , whereas the velocity kernels have units Mm -3 / (km/s). 
Normal modes n = 0 to n— 14 were used to compute the Green’s function and a filter was applied to isolate 
the pi-mode ridge. The cuts are at taken at a height of z — —0.3 Mm. The left panel is a cut along the x axis, 
whereas the right panel is a cut along the y axis. The distance between the two points is 16.0347 Mm. 


where the superscript H is to remind the reader that the assumption of hydrodynamic balance 
has been asserted. Then by equating coefficients we get 

4p( r )=^, P p( r) (43) 

K p,A r ) = K p,M r ) +Po{z)go(z)A(T) (44) 

Equation ( 1431 ) implies that the sound-speed kernel at constant density and pressure in the 
general problem without hydrostatic equilibrium is equal to the sound-speed kernel at con¬ 
stant density in the hydrostatic case. Equation ( 144b demonstrates that more work is required 
to obtain the density kernel at constant sound-speed in the hydrostatic case. It is related to 
the density kernel at constant sound-speed and pressure in the non-hydrodynamic case plus 
an additional term involving the function A(r). Most of the work in going from one set of 
kernels to another involves computing A(r). 


8 Discussion 

We have reviewed the sensitivity of time-distance travel times to small-amplitude steady 
changes in the structure of a solar model as well as their sensitivity to weak and steady 
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flows. In the course of this paper, we introduced a new convenient method for computing 
Green’s functions based on normal-mode summation using the normal modes of the adia¬ 
batic wave equation (in contrast with the approach of lBirch et al .1120041 . based on the modes 
of the damped problem). The sensitivity kernels we computed are essential for time-distance 
inversions of the solar interior. 

It is difficult to disentangle density and pressure perturbations and this will complicate 
time-distance inversions. Our calculations have shown that, to a good approximation, there 
exists a constant a > 0, such that K pcP — —aK Pcp . We have 


St = 




5(2 +K §P + K §P 

— + K-P.cp -y + K P ,cP— 



(45) 


where u = P/p a , K c i u = K c 2 Pp , and K u c = K Pcp . Thus, to a good approximation, two 
thermodynamic quantities are enough to describe the medium. If the density and pressure 
perturbations are in hydrostatic equilibrium, then two thermodynamic quantities are exactly 
enough (see Section[7]). 

The correct treatment of instrumental, geometrical, and mode physics effects is an im¬ 
portant future goal for accurately imaging the solar interior. Correctly including line-of-sight 
and foreshortening effects in the computation of the_kemels, will allow accurate applications 
using large line-of-sight angles, jjackiewicz et al.) 120071) studied line-of-sight effects for 2D 
kernels (for surface waves) and found striking changes to the appearance of kernels that 
will be important for accurately inverting for small-scale flows. Computing travel-time sen¬ 
sitivity kernels in full spherical geometry, rather than using the Cartesian approximation, is 
important for larger distances between the observation points. 

The kernels discussed in this review assume that the Lagrangian pressure perturbation 
vanishes at the top boundary. This is not suitable for frequencies above the acoustic cut-off 
frequency and thus needs to be improved. Another important question is how much do time- 
distance inversion results depend on accurate descriptions of input physics, such as source 
models or damping. 

Most cases assume a laterally homogeneous reference (non-magnetic) model. To study 
more complex features on the Sun (e.g. sunspots), kernels with laterally heterogeneous ref¬ 
erence models are im p ortant. Initial steps tow ards this goal has already been achieved by 
lHanasoge et al.1 J201 lli i lllanasoge et al.1120 1 3) . 
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